library(tidyverse)
library(vdemdata)
library(ggplot2)
library(ggpubr)
library(ggthemes)
library(ggpmisc)
library(GGally)
library(ggrepel)
library(patchwork)
library(extrafont)
library(countrycode)
library(psych)
library(readxl)
library(grid)
library(plm)
library(forcats)
library(stargazer)
library(logr)
theme_set(theme_bw(base_size=11))


## Create a "Data" sub-folder and place all data files in this sub-folder
## Set working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Open log file
log_file <- log_open("replication_log")


###############
## Load data ##
###############
## Main data
vindoc <- readRDS("Data/vindoc_cy.rds") %>%
  mutate(v2xed_ed_dmcon_rec=v2xed_ed_dmcon-0.5,
         v2xed_ed_indcomp=v2xed_ed_inpt*v2xed_ed_dmcon_rec) %>%
  left_join(., vdem %>%
              select(country_id, year, 23:ncol(.))) %>%
  left_join(., read.csv("Data/coder_count.csv", fileEncoding="UTF-8-BOM") )

## Setup map data
world_map=map_data("world") %>% 
  filter(! long > 180) %>% 
  mutate(country_text_id=countrycode(region, origin="country.name", 
                                     destination="iso3c")) %>%
  left_join(., vindoc %>% 
              filter(year==2021) %>% 
              select(country_text_id, 
                     v2xed_ed_inpt,
                     v2xed_ed_inpt_mr,
                     v2xed_ed_dmcon, 
                     v2xed_ed_ptcon,
                     v2xed_ed_con,
                     v2xed_ed_indcomp,
                     v2xedvd_me_inco,
                     coder_count)) %>%
  filter(!(lat<=13 & long<=-115))

## Setup Linz data
linz_dat <- read_xlsx("Data/GWF Autocratic Regimes.xlsx", sheet="TSCS data") %>%
  mutate(gwf_country=str_replace(gwf_country, "Czechoslovakia", "Czech Republic"),
         gwf_country=str_replace(gwf_country, "Cen African Rep", "Central African Republic"),
         gwf_country=str_replace(gwf_country, "Yugoslavia", "Serbia")) %>%
  mutate(country_id=countrycode(gwf_country, origin="country.name", destination="vdem")) %>%
  select(gwf_country, year, country_id, gwf_party, gwf_personal, gwf_monarch, gwf_military) %>%
  left_join(., vindoc %>%
              mutate(lgdppc=log(e_gdppc)))

## Setup coder data
vindoc_coder <- readRDS("Data/vindoc_coder.rds") %>%
  mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  select(country_id, country_name, year, coder_id, contains("conf")) %>%
  mutate(mean_conf=rowMeans(select(., contains("conf")), na.rm=T)) %>%
  left_join(., vdem %>% select(-c(country_name, contains("_conf")))) %>%
  mutate(lgdppc=log(e_gdppc),
         reg_type=ifelse(v2x_regime>=2, 1, 0),
         coup=ifelse(e_pt_coup>0, 1, 0)) %>%
  left_join(., read.csv("Data/is_war_cy.csv") %>%
              mutate(country_name=ifelse(country_name=="Yemen Arab Republic", "Yemen", country_name),
                     country_name=ifelse(country_name=="Czechoslovakia", "Czechia", country_name),
                     country_name=ifelse(country_name=="Yugoslavia", "Serbia", country_name)) %>%
              mutate(country_id=countrycode(country_name, "country.name", "vdem")) %>%
              select(country_id, year, is_war) %>%
              filter(!is.na(country_id))) %>%
  mutate(is_war=ifelse(is.na(is_war) & year<=2007, 0, is_war)) %>%
  left_join(., read.csv("Data/cv_war_cy.csv") %>%
              mutate(country_id=countrycode(CcodeA, "cown", "vdem")) %>%
              mutate(country_id=ifelse(CcodeA==345, 198, country_id),
                     country_id=ifelse(CcodeA==678, 14, country_id)) %>%
              select(-CcodeA)) %>%
  mutate(cv_war=ifelse(cv_war>0, 1, 0),
         cv_war=ifelse(is.na(cv_war) & year<=2014, 0, cv_war)) %>%
  left_join(., read.csv("Data/ert.csv") %>% select(country_id, year, reg_trans))



####################################
## Functions for generating plots ##
####################################

## Generate coefplots for for top/bottom 5 cases
coef_func <- function(plotvar, titletext, ymin, ymax){
  plotvar_high <- paste0(plotvar,"_codehigh68")
  plotvar_low <- paste0(plotvar,"_codelow68")
  plotvar_mr <- paste0(plotvar,"_mr")
  
  if (titletext==""){
    p <- vindoc %>%
      filter(year==2021) %>%
      select(country_name,
             all_of(plotvar),
             all_of(plotvar_high),
             all_of(plotvar_low),
             all_of(plotvar_mr)) %>%
      filter(.[,5]>=3) %>% ## Filtering for number of expert coders
      select(1:4) %>%
      rename("pe"=2,
             "ub"=3,
             "lb"=4) %>%
      arrange(desc(pe)) %>%
      na.omit() %>%
      slice(1:5, (n()-4):n()) %>%
      mutate(country_name=factor(country_name)) %>%
      mutate(country_name=reorder(country_name, pe)) %>%
      ggplot(aes(x=country_name, y=pe))+
      theme(panel.background=element_blank(),
            axis.line=element_line(color="black"),
            plot.title=element_text(size=12, hjust=0.5, color="black"),
            axis.text=element_text(size=12, color="black"),
            axis.text.y=element_text(vjust=0.35),
            axis.title.y=element_blank(),
            axis.title.x=element_blank(),
            legend.title=element_blank(),
            legend.key=element_blank(),
            legend.position="bottom") +
      scale_y_continuous(limits=c(ymin,ymax)) +
      geom_linerange(aes(ymin=lb, ymax=ub), size=1, color="black") +
      geom_point(size=2) +
      geom_vline(aes(xintercept=5.5), size=0.5, color="black", linetype="dashed") +
      coord_flip()
  } else{
    p <- vindoc %>%
      filter(year==2021) %>%
      select(country_name,
             all_of(plotvar),
             all_of(plotvar_high),
             all_of(plotvar_low),
             all_of(plotvar_mr)) %>%
      filter(.[,5]>=3) %>% ## Filtering for number of expert coders
      select(1:4) %>%
      rename("pe"=2,
             "ub"=3,
             "lb"=4) %>%
      arrange(desc(pe)) %>%
      na.omit() %>%
      slice(1:5, (n()-4):n()) %>%
      mutate(country_name=factor(country_name)) %>%
      mutate(country_name=reorder(country_name, pe)) %>%
      ggplot(aes(x=country_name, y=pe))+
      theme(panel.background=element_blank(),
            axis.line=element_line(color="black"),
            plot.title=element_text(size=12, hjust=0.5, color="black"),
            axis.text=element_text(size=12, color="black"),
            axis.text.y=element_text(vjust=0.35),
            axis.title.y=element_blank(),
            axis.title.x=element_blank(),
            legend.title=element_blank(),
            legend.key=element_blank(),
            legend.position="bottom") +
      scale_y_continuous(limits=c(ymin,ymax)) +
      ggtitle(titletext) +
      geom_linerange(aes(ymin=lb, ymax=ub), size=1, color="black") +
      geom_point(size=2) +
      geom_vline(aes(xintercept=5.5), size=0.5, color="black", linetype="dashed") +
      coord_flip()
  }
  
  print(p)
}

## Generate coefplots for all countries
coef_all_func <- function(plotvar){
  plotvar_high <- paste0(plotvar,"_codehigh68")
  plotvar_low <- paste0(plotvar,"_codelow68")
  plotvar_mr <- paste0(plotvar,"_mr")
  
  p_dat <- vindoc %>%
    filter(year==2021) %>%
    select(country_name,
           plotvar,
           plotvar_high,
           plotvar_low,
           plotvar_mr) %>%
    rename("pe"=2,
           "ub"=3,
           "lb"=4,
           "mr"=5) %>%
    arrange(desc(pe)) %>%
    na.omit() %>%
    mutate(country_name=factor(country_name)) %>%
    mutate(country_name=reorder(country_name, pe)) %>%
    mutate(coder=ifelse(mr<3, 0, 1)) %>%
    mutate(coder=factor(coder))
  
  p1 <- p_dat[1:(ceiling(nrow(p_dat)/2)),] %>%
    ggplot(aes(x=country_name, y=pe, color=coder)) +
    theme(panel.background=element_blank(), 
          axis.line=element_line(color="black"),
          plot.title=element_text(size=12, hjust=0.5, color="black"),
          axis.text=element_text(size=11.5, color="black"),
          axis.text.y=element_text(vjust=0.35),
          axis.title.y=element_blank(), 
          axis.title.x=element_blank(), 
          legend.title=element_blank(), 
          legend.key=element_blank(), 
          legend.position="none") +
    scale_y_continuous(limits=c(0,1)) +
    scale_color_manual(values=c("grey60", "black")) +
    geom_linerange(aes(ymin=lb, ymax=ub), size=1) +
    geom_point(size=1.25) +
    coord_flip()
  p2 <- p_dat[(ceiling(nrow(p_dat)/2)+1):nrow(p_dat),] %>%
    ggplot(aes(x=country_name, y=pe, color=coder)) +
    theme(panel.background=element_blank(), 
          axis.line=element_line(color="black"),
          plot.title=element_text(size=12, hjust=0.5, color="black"),
          axis.text=element_text(size=11.5, color="black"),
          axis.text.y=element_text(vjust=0.35),
          axis.title.y=element_blank(), 
          axis.title.x=element_blank(), 
          legend.title=element_blank(), 
          legend.key=element_blank(), 
          legend.position="none") +
    scale_y_continuous(limits=c(0,1)) +
    scale_color_manual(values=c("grey60", "black")) +
    geom_linerange(aes(ymin=lb, ymax=ub), size=1) +
    geom_point(size=1.25) +
    coord_flip()
  p <- ggarrange(p1, p2, nrow=1)
  
  print(p)
}


## Plot variables across regime type
regime_plot <- function(plotvar, titletext, demvar, ymin, ymax, ylimmin, ylimmax, y_break, col_num){
  
  plotvar_high <- paste0(plotvar,"_codehigh68")
  plotvar_low <- paste0(plotvar,"_codelow68")  
  
  if (demvar=="vdem"){
    max_year <-2021
    
    p_dat <- vindoc %>% 
      filter(!is.na(v2x_regime)) %>%
      mutate(regime=ifelse(v2x_regime<2, 0, 1)) %>%
      mutate(regime=case_when(regime==0 ~ "Autocracies", 
                              regime==1 ~ "Democracies")) %>%
      dplyr::select(country_text_id, 
                    year, 
                    regime, 
                    all_of(plotvar),
                    all_of(plotvar_high),
                    all_of( plotvar_low)) %>%
      gather(var, value, -country_text_id, -year, -regime) %>%
      group_by(regime, year, var) %>% 
      summarise(mean.value=mean(value, na.rm=TRUE)) %>%
      mutate(regime=factor(regime, levels=c("Autocracies",
                                            "Democracies")))
  } else{
    max_year <-2015
    
    p_dat <- vindoc %>% 
      filter(!is.na(e_boix_regime)) %>%
      mutate(regime=case_when(e_boix_regime==0 ~ "Autocracies", 
                              e_boix_regime==1 ~ "Democracies")) %>%
      dplyr::select(country_text_id, 
                    year, 
                    regime, 
                    all_of(plotvar),
                    all_of(plotvar_high),
                    all_of( plotvar_low)) %>%
      gather(var, value, -country_text_id, -year, -regime) %>%
      group_by(regime, year, var) %>% 
      summarise(mean.value=mean(value, na.rm=TRUE)) %>%
      mutate(regime=factor(regime, levels=c("Autocracies",
                                            "Democracies")))
  }
   
  
  p_dat$vartype <- titletext[1]
  for (i in 1:nrow(p_dat)){
    for (j in 2:length(plotvar)){
      p_dat$vartype[i] <- ifelse(grepl(plotvar[j], p_dat$var[i], fixed=TRUE), titletext[j], p_dat$vartype[i])
    }
  }
  
  p_dat <- p_dat %>%
    mutate(vartype=factor(vartype)) %>%
    mutate(vartype=fct_relevel(vartype, titletext))
  
  p_dat$var[grep("high", p_dat$var)] <- "ub"
  p_dat$var[grep("low", p_dat$var)] <- "lb"
  p_dat$var[!grepl("lb", p_dat$var) & !grepl("ub", p_dat$var)] <- "pe"
  
  p_dat <- p_dat %>%
    pivot_wider(names_from=var,
                values_from=mean.value)
  
  p <- p_dat %>% 
    ggplot(aes(x=year, y=pe, fill=regime, color=regime)) + 
    geom_line(aes(linetype=regime), lwd=1) +
    geom_ribbon(aes(ymin=lb, ymax=ub), alpha=0.5, color=NA) +
    scale_x_continuous(breaks=seq(1945,2015,10), limits=c(1945,max_year), expand=c(0,0)) +
    scale_y_continuous(breaks=seq(ymin,ymax,y_break), limits=c(ylimmin,ylimmax)) +
    scale_linetype_manual(values=c("dotted", "solid")) +
    guides(fill=guide_legend(ncol=2),
           color=guide_legend(override.aes=list(fill=NA))) +
    theme(legend.position="bottom", 
          legend.title=element_blank(),
          axis.text.x=element_text(size=rel(1.5)),
          axis.text.y=element_text(size=rel(1.5)),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          legend.text=element_text(size=rel(1.5)),
          strip.background=element_rect(fill="white"),
          strip.text.x=element_text(size=rel(1.5), color="black"))+
    facet_wrap(~vartype, ncol=col_num)
  
  p <- p + theme(panel.spacing=unit(1.5, "lines"))
  
  print(p)
}


## Generate functions used in paired plots
cor_func <- function(data, mapping, method, ...){
  x <- eval_data_col(data, mapping$x)
  y <- eval_data_col(data, mapping$y)
  
  corr <- cor(x, y, method=method, use="pairwise.complete.obs")
  
  ggally_text(
    label=as.character(round(corr, 2)),
    mapping=aes(),
    xP=0.5, yP=0.5,
    color='black',
    ...
  )
}

loess_pt_fn <- function(data, mapping, method="loess", ...){
  p <- ggplot(data=data, mapping=mapping) +
    geom_point() +
    geom_smooth(method=method, ...)
  p
}


## Plot for generating correlations between democracy measures and V-Indoc indices
corr_dem <- function(vindoc_var, vindoc_var_lab){
  corr_dat <- vindoc %>%
    mutate(e_p_polity=ifelse(-10>e_p_polity, 2, e_p_polity)) %>%
    select(country_id, year, v2x_liberal, v2x_polyarchy, e_p_polity, e_uds_median, all_of(vindoc_var)) %>%
    filter(year==2010) %>%
    select(-c(country_id, year)) %>%
    rename("Liberal Democracy"=1,
           "Electoral Democracy"=2,
           "Polity 5"=3,
           "UDS"=4,
           "vindoc"=5) %>%
    gather(key="variable", value="value", 1:4) %>%
    mutate(variable=factor(variable, levels=c("Liberal Democracy",
                                              "Electoral Democracy",
                                              "Polity 5",
                                              "UDS")))
  
  corr_out <- ggplot(corr_dat, aes(x=value, y=vindoc)) +
    theme(panel.background=element_blank(), 
          axis.line=element_line(color="black"),
          plot.title=element_text(size=12, hjust=0.5, color="black"),
          axis.text=element_text(size=12, color="black"),
          axis.text.y=element_text(size=12, color="black"),
          axis.title.y=element_text(size=12, color="black"), 
          axis.title.x=element_blank(), 
          legend.title=element_blank(), 
          legend.key=element_blank(), 
          legend.position="none",
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          panel.spacing=unit(0.2, "cm")) +
    geom_point() +
    geom_smooth() +
    ylab(vindoc_var_lab) +
    stat_cor(aes(label=..r.label..)) +
    scale_y_continuous(limits=c(0,1.1), breaks=c(0,1)) +
    facet_wrap(~variable, 
               nrow=1,
               scales="free_x")
  
  return(corr_out)
}


## Linz analysis and plot
linz_out <- function(var, coder_ct_var, model){
  
  print_order <- c(var,
                   "gwf_party",
                   "gwf_personal",
                   "gwf_military",
                   "gwf_monarch",
                   "lgdppc")
  
  linz_subset <- linz_dat %>%
    select(country_id,
           year,
           all_of(print_order),
           all_of(coder_ct_var)) %>%
    mutate(sample_full=ifelse(rowSums(is.na(.))==0, 1, 0),
           sample_3coder=ifelse(.[[coder_ct_var]]>=3, 1, 0)) %>%
    select(-all_of(coder_ct_var))
  
  ## Summary statistics
  sum_full <- linz_subset %>%
    filter(sample_full==1) %>%
    select(all_of(print_order)) %>%
    pivot_longer(
      everything(),
      names_to="Variable",
      values_to="Value"
    ) %>%
    mutate(Variable=factor(Variable, levels=print_order)) %>%
    group_by(Variable) %>%
    summarize(
      N=n(),
      Mean=round(mean(Value, na.rm=T), 4),
      SD=round(sd(Value, na.rm=T), 4),
      Min=round(min(Value, na.rm=T), 4),
      Max=round(max(Value, na.rm=T), 4)) %>%
    as.data.frame()
  
  sum_3coder <- linz_subset %>%
    filter(sample_3coder==1) %>%
    select(all_of(print_order)) %>%
    pivot_longer(
      everything(),
      names_to="Variable",
      values_to="Value"
    ) %>%
    mutate(Variable=factor(Variable, levels=print_order)) %>%
    group_by(Variable) %>%
    summarize(
      N=n(),
      Mean=round(mean(Value, na.rm=T), 4),
      SD=round(sd(Value, na.rm=T), 4),
      Min=round(min(Value, na.rm=T), 4),
      Max=round(max(Value, na.rm=T), 4)) %>%
    as.data.frame()
  
  ## Run FE regressions
  fe_reg <- linz_subset %>%
    filter(sample_full==1) %>%
    plm(as.formula(paste0(var," ~ gwf_party + gwf_personal + gwf_monarch + lgdppc")),
        data=., 
        index=c("country_id", "year"), 
        model="within",
        effect=model)
  
  fe_reg_3coder <- linz_subset %>%
    filter(sample_3coder==1) %>%
    plm(as.formula(paste0(var," ~ gwf_party + gwf_personal + gwf_monarch + lgdppc")), 
        data=., 
        index=c("country_id", "year"), 
        model="within",
        effect=model)
  
  ## Coefplots
  plot_full <- as.data.frame(cbind(fe_reg$coefficients, confint(fe_reg))) %>%
    rename(coef=1,
           lb=2,
           ub=3) %>%
    mutate(vars=factor(c("Party", "Personalist", "Monarchy", "GDPpc (log)"),
                       levels=c("GDPpc (log)", "Monarchy", "Personalist",  "Party"))) %>%
    mutate(subset="All observations")
  
  plot_3coder <- as.data.frame(cbind(fe_reg_3coder$coefficients, confint(fe_reg_3coder))) %>%
    rename(coef=1,
           lb=2,
           ub=3) %>%
    mutate(vars=factor(c("Party", "Personalist", "Monarchy", "GDPpc (log)"),
                       levels=c("GDPpc (log)", "Monarchy", "Personalist",  "Party"))) %>%
    mutate(subset="Observations with at least 3 coders")
  
  plot_dat <- rbind(plot_full, plot_3coder) %>%
    mutate(subset=factor(subset))
  
  cp <- ggplot(plot_dat, aes(x=vars, y=coef, group=fct_rev(subset), color=subset)) +
    theme(axis.title.y=element_blank(),
          axis.text.y=element_text(size=10, color="black"),
          axis.ticks.y=element_blank(),
          axis.title.x=element_blank(),
          axis.text.x=element_text(size=10, color="black"),
          legend.position="bottom",
          legend.title=element_blank(),
          legend.key=element_rect(fill=NA),
          legend.text=element_text(size=10),
          panel.border=element_rect(colour="black", fill=NA, linewidth=1),
          panel.background=element_rect(fill="transparent"),
          panel.grid.major=element_line(color="lightgrey"),
          panel.grid.minor=element_line(color="lightgrey")) +
    geom_errorbar(aes(ymin=lb, ymax=ub, width=0), position=position_dodge(width=0.45), size=1) +
    geom_point(aes(shape=subset), position=position_dodge(width=0.45), size=2) +
    geom_hline(yintercept=0, linetype="dashed") +
    guides(shape=guide_legend(ncol=1)) +
    scale_y_continuous(limits=c(-0.2,0.1)) +
    scale_shape_manual(values=c(19, 17)) +
    coord_flip()
  
  out_list <- list(sum_full=sum_full,
                   sum_3coder=sum_3coder,
                   fe_reg=fe_reg,
                   fe_reg_3coder=fe_reg_3coder,
                   plot_dat=plot_dat,
                   cp=cp)
  
  return(out_list)
}

## Plot function for coder confidence across variables
coder_event_out <- function(event_var, is_dummy, dummy_level, var_lab){
  if (is_dummy==0){
    plot_dat <- vindoc_coder %>%
      select(year, all_of(event_var), mean_conf) %>%
      na.omit %>%
      rename(var=2)
    
    plot_out <- ggplot(plot_dat, aes(x=var, y=mean_conf)) +
      theme(plot.title=element_text(size=12, hjust=0.5),
            axis.title.y=element_blank(),
            axis.text.y=element_text(size=12, color="black"),
            axis.title.x=element_blank(),
            axis.text.x=element_text(size=12, color="black"),
            legend.title=element_blank(),
            legend.position="none",
            legend.key=element_rect(fill=NA),
            legend.text=element_text(size=12),
            panel.border=element_rect(colour = "black", fill=NA, size=1),
            panel.background=element_rect(fill="white"),
            panel.grid=element_line(color="grey"),
            panel.grid.minor=element_blank(),
            plot.caption=element_text(hjust=0, size=12))+
      geom_point(shape=16, alpha=0.2, size=0.5) +
      # geom_smooth() +
      ggtitle(var_lab) +
      scale_y_continuous(limits=c(0, 1), breaks=round(seq(0,1,0.2),1))
    
  } else{
    plot_dat <- vindoc_coder %>%
      select(year, all_of(event_var), mean_conf) %>%
      na.omit %>%
      rename(var=2) %>%
      mutate(var=factor(var, labels=dummy_level))
    
    
    plot_out <- ggplot(plot_dat, aes(x=var, y=mean_conf)) +
      theme(plot.title=element_text(size=12, hjust=0.5),
            axis.title.y=element_blank(),
            axis.text.y=element_text(size=12, color="black"),
            axis.title.x=element_blank(),
            axis.text.x=element_text(size=10, color="black"),
            legend.title=element_blank(),
            legend.position="none",
            legend.key=element_rect(fill=NA),
            legend.text=element_text(size=12),
            panel.border=element_rect(colour = "black", fill=NA, size=1),
            panel.background=element_rect(fill="white"),
            panel.grid=element_line(color="grey"),
            panel.grid.minor=element_blank(),
            plot.caption=element_text(hjust=0, size=12))+
      geom_boxplot(outlier.size=0.5) +
      ggtitle(var_lab) +
      scale_y_continuous(limits=c(0, 1), breaks=round(seq(0,1,0.2),1))
  }
  
  return(plot_out)
}


## Indices across democracies
dem_splot <- function(var, var_lab){
  plot_dat <- vindoc %>%
    select(country_text_id, year, v2x_libdem, all_of(var)) %>%
    filter(year==2021) %>%
    select(-year) %>%
    rename(libdem=v2x_libdem,
           vindoc_var=3)
   
  p1 <- ggplot(plot_dat, aes(x=libdem, y=vindoc_var, label=country_text_id)) +
    geom_point() +
    geom_text_repel(box.padding=0.2, point.padding=0.1, size=4, max.overlaps=30, color="black") +
    scale_x_continuous(limits=c(0,1), breaks=c(0,0.5,1)) +
    scale_y_continuous(limits=c(-0.075,1.1), breaks=c(0,0.5,1)) +
    theme_bw() +
    theme(axis.title.x=element_text(size=14),
          axis.title.y=element_text(size=14),
          axis.text.x=element_text(size=14, color="black"),
          axis.text.y=element_text(size=14, color="black")) +
    xlab("Liberal Democracy Index") +
    ylab(var_lab) +
    annotate("text", x=0.0475, y=1.075, label=paste0("Corr=", round(cor(plot_dat[,2], plot_dat[,3], use="pairwise.complete.obs"),2)), size=6)
    
    return(p1)
}



##################################################
## Figure 4: Number of unique coders by country ##
##################################################
fig4 <- world_map %>%
  mutate(coder_count=ifelse(is.na(coder_count), -1, coder_count)) %>%
  mutate(coder_count=cut(coder_count, breaks=c(-Inf,0,3,5,10,Inf), labels=c("N/A","[1,3)","[3,5)","[5,10)","[10,20]"))) %>%
  mutate(coder_count=factor(coder_count)) %>%
  # ggplot(aes(fill=coder_count, map_id=region, alpha=is.na(coder_count))) +
  ggplot(aes(fill=coder_count, map_id=region)) +
  geom_map(map=world_map, col=NA) +
  coord_fixed(1.3) +
  expand_limits(x=world_map$long, y=world_map$lat) +
  ggthemes::theme_map() +
  # scale_fill_brewer(breaks=c("N/A", "[1,3)","[3,5)","[5,10)","[10,20]"),
  #                   na.value="gray60", palette="Reds") +
  # scale_alpha_manual(values = c("TRUE" = .5), labels = c("N/A")) +
  scale_fill_manual(values=c("grey60", "#FEE5D9", "#FCAE91", "#FB6A4A", "#CB181D")) +
  # scale_alpha_manual(values=c(0.5,1,1,1,1), guide=F) +
  scale_y_continuous(limits=c(-50,90)) +
  theme(legend.text=element_text(size=rel(2.5)),
        legend.title=element_blank(),
        legend.position=c(0.05,0.015),
        plot.title=element_text(size=rel(2), hjust=0.5),
        plot.background=element_rect(color="black", fill="white"),
        plot.caption=element_text(hjust=0, size=rel(1.45))) +
  guides(fill=guide_legend(reverse=T))

ggsave(fig4, file="Figures/Figure 4.pdf", height=6, width=12)



##################################################################################
## Figure 5: Percentage and number of countries covered in the V-Indoc data set ##
##################################################################################
total_countries <- vdem %>% distinct(country_text_id, year) %>%
  group_by(year) %>% summarise(total=n()) %>% filter(year>=1945)

## Sample education question (the first one)

coverage.5coders.educ <- vindoc %>% filter(v2edcentcurrlm_nr>=5) %>%
  # mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  distinct(country_text_id, year) %>% ungroup() %>%
  left_join(., total_countries) %>%
  group_by(year) %>%
  summarise(share=n()/total, n=n()) %>% mutate(coders="at least five coders")
coverage.3coders.educ <- vindoc %>% filter(v2edcentcurrlm_nr>=3) %>%
  # mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  distinct(country_text_id, year) %>% ungroup() %>%
  left_join(., total_countries) %>%
  group_by(year) %>%
  summarise(share=n()/total, n=n()) %>% mutate(coders="at least three coders")
coverage.1coder.educ <- vindoc %>% filter(v2edcentcurrlm_nr>=1) %>%
  # mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  distinct(country_text_id, year) %>% ungroup() %>%
  left_join(., total_countries) %>%
  group_by(year) %>%
  summarise(share=n()/total, n=n()) %>% mutate(coders="at least one coder")

coders.educ <-
  coverage.5coders.educ %>% bind_rows(., coverage.3coders.educ) %>%
  bind_rows(., coverage.1coder.educ)

coders.educ$coders <- factor(coders.educ$coders, levels = c('at least five coders',
                                                            'at least three coders',
                                                            'at least one coder'))

## Sample media question (the first one)

coverage.5coders.media <- vindoc %>% filter(v2medstateprint_nr>=5) %>%
  # mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  distinct(country_text_id, year) %>% ungroup() %>%
  left_join(., total_countries) %>%
  group_by(year) %>%
  summarise(share=n()/total, n=n()) %>% mutate(coders="at least five coders")
coverage.3coders.media <- vindoc %>% filter(v2medstateprint_nr>=3) %>%
  # mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  distinct(country_text_id, year) %>% ungroup() %>%
  left_join(., total_countries) %>%
  group_by(year) %>%
  summarise(share=n()/total, n=n()) %>% mutate(coders="at least three coders")
coverage.1coder.media <- vindoc %>% filter(v2medstateprint_nr>=1) %>%
  # mutate(year=as.numeric(substr(historical_date, 1, 4))) %>%
  distinct(country_text_id, year) %>% ungroup() %>%
  left_join(., total_countries) %>%
  group_by(year) %>%
  summarise(share=n()/total, n=n()) %>% mutate(coders="at least one coder")

coders.media <-
  coverage.5coders.media %>% bind_rows(., coverage.3coders.media) %>%
  bind_rows(., coverage.1coder.media)

coders.media$coders <- factor(coders.media$coders, levels = c('at least five coders',
                                                              'at least three coders',
                                                              'at least one coder'))

coders <- bind_rows(coders.educ %>% mutate(label = "Education"),
                    coders.media %>% mutate(label = "Media")) %>%
  gather(var, value, -year, -coders, -label)



fig5 <- coders %>%
  mutate(coders=factor(coders, levels=c("at least five coders",
                                        "at least three coders",
                                        "at least one coder"))) %>%
  ggplot() + 
  geom_line(data = subset(coders, var=="share"), aes(x=year, y=value, color=coders, size=coders)) + 
  geom_line(data = subset(coders, var=="n"), 
            aes(x=year, y=value/160, color=coders, size=coders), lty="dotted") +
  facet_wrap(~ label, scales = "free_y") +
  scale_y_continuous("Percentage of countries", limits=c(.2,1), 
                     labels = scales::percent,
                     sec.axis = sec_axis(trans = (~(.-0)*160), name = "Number of countries", 
                                         breaks=seq(0, 170, 20))
  ) +
  scale_size_manual(values=c(1.4, 1, 0.6)) +
  scale_x_continuous(breaks=seq(1945,2021,10), limits=c(1945,2021)) +
  labs(x="Year") + 
  guides(fill=guide_legend(ncol=2)) +
  theme(legend.position="bottom", legend.title=element_blank(),
        axis.text.x = element_text(size=rel(1), angle = 90, vjust = 0.5, hjust=1),
        axis.text.y=element_text(size=rel(1)),
        axis.title.x=element_text(size=rel(1)),
        axis.title.y=element_text(size=rel(1)),
        legend.text=element_text(size=rel(1)),
        strip.background = element_rect(fill="white"),
        strip.text.x = element_text(
          size = 12, color = "black"),
        axis.title=element_blank()) +
  # scale_color_viridis("", discrete=T) +
  # scale_fill_viridis("", discrete=T) +
  guides(color=guide_legend(nrow=1, byrow=TRUE), alpha = FALSE) 

ggsave(fig5, file="Figures/Figure 5.pdf", height=4, width=7)



############################################################
## Figure 6: Indoctrination potential in education (2021) ##
############################################################
fig6 <- world_map %>%
  ggplot(aes(fill=v2xed_ed_inpt, map_id=region, alpha = is.na(v2xed_ed_inpt))) +
  geom_map(map=world_map, col=NA) +
  coord_fixed(1.3) +
  expand_limits(x=world_map$long, y=world_map$lat) +
  ggthemes::theme_map() +
  scale_fill_continuous("Potential: ", 
                        low="thistle1",
                        high="darkred",
                        guide="colorbar", 
                        na.value="grey60",
                        limits=c(0,1),
                        breaks=c(0,0.5,1)
  ) +
  scale_alpha_manual(values = c("TRUE" = 0.5), labels = c("N/A")) +
  scale_y_continuous(limits=c(-50,90)) +
  theme(legend.text=element_text(size=rel(2.5)),
        legend.title=element_blank(),
        legend.position=c(0.05,0.015),
        plot.title=element_text(size=rel(2), hjust=0.5),
        plot.background=element_rect(color="black", fill="white"),
        plot.caption=element_text(hjust=0, size=rel(1.45)))

ggsave(fig6, file="Figures/Figure 6.pdf", height=6, width=12)




##############################################################################
## Figure 7: Indoctrination potential in education (2021): Bottom/top cases ##
##############################################################################
fig7 <- coef_func("v2xed_ed_inpt", "", 0, 1)
ggsave(fig7, file="Figures/Figure 7.pdf", height=4, width=5.5)



##########################################################
## Figure 8: Indoctrination content in education (2021) ##
##########################################################
fig8_1 <- world_map %>% 
  ggplot(aes(fill=v2xed_ed_dmcon, map_id=region, alpha = is.na(v2xed_ed_dmcon))) +
  geom_map(map=world_map, col=NA) +
  coord_fixed(1.3) +
  expand_limits(x=world_map$long, y=world_map$lat) +
  ggthemes::theme_map() +
  scale_fill_continuous("Content: ", 
                        low="white",
                        high="dodgerblue4",
                        guide="colorbar", 
                        na.value="gray60",
                        limits=c(0,1),
                        breaks=c(0,0.5,1)
  ) +
  scale_alpha_manual(values = c("TRUE" = .5), labels = c("N/A")) +
  scale_y_continuous(limits=c(-50,90)) +
  ggtitle("Democratic Content") +
  theme(legend.text=element_text(size=rel(2)),
        legend.title=element_blank(),
        legend.position="none",
        plot.title=element_text(size=rel(2.5), hjust=0.5),
        plot.background=element_rect(color="black", fill="white"),
        plot.caption=element_text(hjust=0, size=rel(1.45)))

fig8_2 <- world_map %>% 
  ggplot(aes(fill=v2xed_ed_ptcon, map_id=region, alpha = is.na(v2xed_ed_ptcon))) +
  geom_map(map=world_map, col=NA) +
  coord_fixed(1.3) +
  expand_limits(x=world_map$long, y=world_map$lat) +
  ggthemes::theme_map() +
  scale_fill_continuous("Content: ", 
                        low="white",
                        high="dodgerblue4",
                        guide="colorbar", 
                        na.value="gray60",
                        limits=c(0,1),
                        breaks=c(0,0.5,1)
  ) +
  scale_alpha_manual(values = c("TRUE" = .5), labels = c("N/A")) +
  scale_y_continuous(limits=c(-50,90)) +
  ggtitle("Patriotic Content") +
  theme(legend.text=element_text(size=rel(2)),
        legend.title=element_blank(),
        legend.position=c(0.05,0.015),
        plot.title=element_text(size=rel(2.5), hjust=0.5),
        plot.background=element_rect(color="black", fill="white"),
        plot.caption=element_text(hjust=0, size=rel(1.45)))

fig8 <- ggarrange(fig8_1, fig8_2, nrow=2, ncol=1)
ggsave(fig8, file="Figures/Figure 8.pdf", height=12, width=12)



############################################################################
## Figure 9: Indoctrination content in education (2021): Bottom/top cases ##
############################################################################
fig9_1 <- coef_func("v2xed_ed_dmcon", "Democratic Content Index", 0, 1)
fig9_2 <- coef_func("v2xed_ed_ptcon", "Patriotic Content Index", 0, 1)

fig9 <- ggarrange(fig9_1, fig9_2, nrow=1, ncol=2, widths=c(1.12,1))
ggsave(fig9, file="Figures/Figure 9.pdf", height=4, width=8)



#################################################################
## Figure 10: Democracy and the indoctrination indices in 2021 ##
#################################################################
fig10 <- vindoc %>%
  select(country_id, year, v2x_libdem, v2xed_ed_inpt, v2xed_ed_dmcon, v2xed_ed_ptcon) %>%
  filter(year==2015) %>%
  select(-c(country_id, year)) %>%
  rename("Liberal Democracy"=1, 
         "Indoctrination Potential"=2,
         "Democratic Content"=3
         , "Patriotic Content"=4) %>%
  ggpairs(.,
          diag=list(continuous=wrap("barDiag", color="black")),
          lower=list(continuous=loess_pt_fn),
          upper=list(continuous=wrap(cor_func, size=8, method="pearson"))) +
  theme_bw() +
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())

fig10[1,1] <- fig10[1,1] + 
  ylim(0, 15)
fig10[2,2] <- fig10[2,2] + 
  ylim(0, 15)
fig10[3,3] <- fig10[3,3] +
  ylim(0, 15)
fig10[4,4] <- fig10[4,4] + 
  ylim(0, 15)

fig10[2,1] <- fig10[2,1] + 
  scale_y_continuous(limits=c(0,1), breaks=c(0,1))
fig10[3,1] <- fig10[3,1] + 
  scale_y_continuous(limits=c(0,1), breaks=c(0,1))
fig10[4,1] <- fig10[4,1] + scale_x_continuous(limits=c(0,1), breaks=c(0,1)) + 
  scale_y_continuous(limits=c(0,1), breaks=c(0,1))


fig10[4,2] <- fig10[4,2] + 
  scale_x_continuous(limits=c(0,1), breaks=c(0,1))
fig10[4,3] <- fig10[4,3] + 
  scale_x_continuous(limits=c(0,1), breaks=c(0,1))
fig10[4,4] <- fig10[4,4] + 
  scale_x_continuous(limits=c(0,1), breaks=c(0,1))

ggsave(fig10, file="Figures/Figure 10.pdf", height=8, width=8)



#################################################################################
## Figure 11: Indoctrination potential and content in education across regimes ##
#################################################################################
fig11 <- regime_plot(c("v2xed_ed_inpt",
                   "v2xed_ed_dmcon", 
                   "v2xed_ed_ptcon"),
                 c("Indoctrination potential",
                   "Democratic indoctrination content",
                   "Patriotic indoctrination content"),
                 "vdem",
                 0,1,0,1,0.25,3)

ggsave("Figures/Figure 11.pdf", height=4, width=12)



###########################################################################
## Figure 12: Indoctrination potential and content in education (Russia) ##
###########################################################################
werd <- read_xls("Data/werd.xls") %>%
  mutate(country_text_id = countrycode(countryname, origin = "country.name", 
                                       destination = "iso3c"))

russia.reforms <- werd %>% filter(country_text_id=="RUS") %>% rename(year_reform=year)
tajikistan.reforms <- werd %>% filter(country_text_id=="TJK") %>% rename(year_reform=year)

## Indicators to illustrate aggregate indices

pind_dat1 <-
  vindoc %>%
  select(country_name, 
         year,
         v2edtefire, 
         v2edtefire_codehigh68, 
         v2edtefire_codelow68) %>%
  mutate(var="Political Teacher Firing") %>%
  filter(country_name %in% c("Russia")) %>%
  rename(pe=v2edtefire,
         ub=v2edtefire_codehigh68,
         lb=v2edtefire_codelow68) 

pind_dat2 <-
  vindoc %>%
  select(country_name, 
         year,
         v2edcritical, 
         v2edcritical_codehigh68, 
         v2edcritical_codelow68) %>%
  mutate(var="Critical Discussion in Classroom") %>%
  filter(country_name%in% c("Russia")) %>%
  rename(pe=v2edcritical,
         ub=v2edcritical_codehigh68,
         lb=v2edcritical_codelow68)

pind_dat3 <-
  vindoc %>%
  select(country_name, 
         year,
         v2edscpatriotcb, 
         v2edscpatriotcb_codehigh68, 
         v2edscpatriotcb_codelow68) %>%
  mutate(var="Celebration of Patriotic Symbols") %>%
  filter(country_name %in% c("Russia")) %>%
  rename(pe=v2edscpatriotcb,
         ub=v2edscpatriotcb_codehigh68,
         lb=v2edscpatriotcb_codelow68)

pind_dat <- rbind(pind_dat1, pind_dat2, pind_dat3) %>%
  mutate(var=factor(var)) %>%
  mutate(var=fct_relevel(var, c("Political Teacher Firing", "Critical Discussion in Classroom", "Celebration of Patriotic Symbols")))

## Aggregate indices
p_dat1 <- vindoc %>%
  select(country_name, 
         year,
         v2xed_ed_inpt, 
         v2xed_ed_inpt_codehigh68, 
         v2xed_ed_inpt_codelow68) %>%
  mutate(var="Indoctrination Potential") %>%
  filter(country_name %in% c("Russia")) %>%
  rename(pe=v2xed_ed_inpt,
         ub=v2xed_ed_inpt_codehigh68,
         lb=v2xed_ed_inpt_codelow68)

p_dat2 <- vindoc %>%
  select(country_name, 
         year,
         v2xed_ed_dmcon, 
         v2xed_ed_dmcon_codehigh68, 
         v2xed_ed_dmcon_codelow68) %>%
  mutate(var="Democratic Content") %>%
  filter(country_name %in% c("Russia")) %>%
  rename(pe=v2xed_ed_dmcon,
         ub=v2xed_ed_dmcon_codehigh68,
         lb=v2xed_ed_dmcon_codelow68)

p_dat3 <- vindoc %>%
  select(country_name, 
         year,
         v2xed_ed_ptcon, 
         v2xed_ed_ptcon_codehigh68, 
         v2xed_ed_ptcon_codelow68) %>%
  mutate(var="Patriotic Content") %>%
  filter(country_name %in% c("Russia")) %>%
  rename(pe=v2xed_ed_ptcon,
         ub=v2xed_ed_ptcon_codehigh68,
         lb=v2xed_ed_ptcon_codelow68)

p_dat <- rbind(p_dat1, p_dat2, p_dat3) %>%
  mutate(var=factor(var)) %>%
  mutate(var=fct_relevel(var, c("Indoctrination Potential", "Democratic Content", "Patriotic Content")))

fig12 <-
  p_dat %>% filter(country_name=="Russia") %>%
  ggplot(aes(x=year, y=pe)) + 
  geom_line(lwd=1) +
  geom_ribbon(aes(ymin=lb, ymax=ub), alpha=0.3, color="grey") +
  scale_x_continuous(breaks=seq(1945,2015,10), limits=c(1945,2021), expand=c(0,0)) +
  scale_y_continuous(breaks=seq(0,1,1), limits=c(0,1)) +
  guides(fill=guide_legend(ncol=1),
         color=guide_legend(override.aes=list(fill=NA))) +
  geom_vline(xintercept = russia.reforms$year_reform, color="red", lty=2, alpha=.5) +
  theme(legend.position="none", legend.title=element_blank(),
        axis.text.x=element_text(angle = 90, size=rel(1.5)),
        axis.text.y=element_text(size=rel(1.5)),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        plot.title = element_text(size=rel(1.5)),
        legend.text=element_text(size=rel(1.5)),
        strip.background=element_rect(fill="white"),
        strip.text.x=element_text(size=rel(1.5), color="black"))+
  facet_wrap(~var, ncol=3) +
  
  pind_dat %>% filter(country_name=="Russia") %>%
  ggplot(aes(x=year, y=pe)) + 
  geom_line(lwd=1) +
  geom_ribbon(aes(ymin=lb, ymax=ub), alpha=0.3, color="grey") +
  scale_x_continuous(breaks=seq(1945,2015,10), limits=c(1945,2021), expand=c(0,0)) +
  scale_y_continuous(breaks=seq(-3,3,1), limits=c(-3,3)) +
  guides(fill=guide_legend(ncol=1),
         color=guide_legend(override.aes=list(fill=NA))) +
  geom_vline(xintercept = russia.reforms$year_reform, color="red", lty=2, alpha=.5) +
  theme(legend.position="none", legend.title=element_blank(),
        axis.text.x=element_text(angle = 90, size=rel(1.5)),
        axis.text.y=element_text(size=rel(1.5)),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        plot.title = element_text(size=rel(1.5)),
        legend.text=element_text(size=rel(1.5)),
        strip.background=element_rect(fill="white"),
        strip.text.x=element_text(size=rel(1.5), color="black"))+
  facet_wrap(~var, ncol=3) +
  
  plot_layout(ncol=1)

ggsave(fig12, file="Figures/Figure 12.pdf", height=8, width=12)



#####################################################################################
## Figure 13: Indoctrination potential in education across autocratic regime types ##
#####################################################################################
fig13 <- linz_out(
  "v2xed_ed_inpt",
  "v2xed_ed_inpt_mr",
  "individual")

ggsave(fig13$cp, file="Figures/Figure 13.pdf", height=3, width=6)



################################################################################################################################################
## Figure B-1: Number and proportion of journal articles, books and book chapters that mention the term "indoctrinat*" at least once on JStor ##
################################################################################################################################################
## Data obtained from https://tdm-pilot.org/
dta <- cbind(seq(1900, 2010, 10),
             #economics
             c(1, 1, 0, 3, 30, 60, 80, 114,
               104, 105, 200, 196),
             # education
             c(1, 2, 8, 80, 158, 261, 309,
               294, 236, 308, 374, 332),
             # history
             c(5, 10, 7, 18, 28, 70, 94, 99, 117,
               159, 271, 315),
             # philosophy
             c(19, 33, 54, 329, 793, 1307, 1873,
               2356, 2610, 3720, 5161, 6061),
             # polisci
             c(5, 7, 21, 46, 193, 249, 307,
               349, 362, 522, 1045, 1426),
             # socialsci
             c(5, 6, 4, 40, 109, 218, 270, 379,
               337, 544, 810, 825),
             # total number of texts with indoctrinat*
             c(32, 51, 84, 462, 1189, 1952, 2641, 3199, 3408,
               4781, 6963, 8127),
             # total academic texts on JSTOR
             c(22648, 33556, 39315, 52596,
               68544, 97224, 144976, 229568,
               291867, 380183, 502355, 510325)
) %>% as.data.frame() %>%
  rename(year = V1, n.econ = V2, n.educ = V3, n.hist = V4,
         n.phil = V5, n.polsci = V6, n.socsci = V7,
         n.total.indoctrinat = V8, n.jstor = V9)

dta <- dta %>% mutate(total.check=n.econ+n.educ+n.hist+
                        n.phil+n.polsci+n.socsci,
                      delta = total.check-n.total.indoctrinat)

## Main figure with number and proportion of texts
dta %>%
  ggplot(aes(x=year, y=n.total.indoctrinat, fill="red4",
             color="red4")) +
  geom_col(position = "dodge2") + 
  scale_x_continuous(breaks=seq(1900, 2010, 10)) +
  scale_y_continuous(breaks=seq(0, 9000, 1000)) +
  theme(legend.position="none",
        panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  scale_colour_manual("", values=c("red4")) +
  scale_fill_manual("", values=c("red4")) +
  labs(x="decade", title="(a) Number of works",
       y="Number of texts that mention indoctrinat*") +
  dta %>% 
  mutate(prop.indoctr = n.total.indoctrinat/n.jstor) %>%
  ggplot(aes(x=year, y=prop.indoctr, fill="red4",
             color="red4")) +
  geom_col(position = "dodge2") + 
  scale_x_continuous(breaks=seq(1900, 2010, 10)) +
  theme(legend.position="none",
        panel.background=element_blank(),
        panel.border=element_blank(),
        panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  scale_colour_manual("", values=c("red4")) +
  scale_fill_manual("", values=c("red4")) +
  scale_y_continuous(labels=scales::percent) +
  labs(x="decade", title="(b) Proportion of all publushed works",
       y="% of texts that mention indoctrinat*")

ggsave(file="Figures/Figure B-1.pdf", width=9, height=5)



######################################################################
## Figures B-2 through B-4 generated through https://tdm-pilot.org/ ##
######################################################################



####################################################################################
## Figure B-5: Levels of democracy and critical engagement with education Content ##
####################################################################################
figb_5 <- ggplot(vindoc, aes(x=v2x_libdem, y=v2edcritical)) +
  theme(plot.title=element_text(size=14),
        axis.title.y=element_text(size=12),
        axis.text.y=element_text(size=12, color="black"),
        axis.title.x=element_text(size=12),
        axis.text.x=element_text(size=12, color="black"),
        legend.title=element_blank(),
        legend.key=element_blank(),
        legend.text=element_text(size=15),
        legend.position="none",
        panel.border=element_rect(colour = "black", fill=NA, size=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey"),
        panel.grid.minor=element_blank(),
        plot.margin=unit(c(0.25,0.25,0.25,0.25),"cm")) +
  geom_point(stroke=NA) +
  stat_cor(aes(label=..r.label..)) +
  xlab("Liberal Democracy Index") +
  ylab("Critical Engagement with Education Content") +
  scale_y_continuous(limits=c(-3,4.35), breaks=seq(-3,4,1)) +
  scale_x_continuous(limits=c(0,1), breaks=c(0, 0.5, 1)) +
  geom_smooth()

ggsave(figb_5, file="Figures/Figure B-5.pdf", height=5, width=5)



##########################################################
## Figure F-9: Centralized curriculum (model estimates) ##
##########################################################
p_dat <- vindoc %>%
  select(country_name, year, v2edcentcurrlm, v2edcentcurrlm_codehigh68, v2edcentcurrlm_codelow68) %>%
  filter(country_name=="Bolivia" | country_name=="South Africa" | country_name=="United States of America")

figf_9 <- p_dat %>% 
  ggplot(aes(x=year, y=v2edcentcurrlm, fill=country_name, color=country_name)) + 
  geom_line(lwd=1) +
  geom_ribbon(aes(ymin=v2edcentcurrlm_codelow68, ymax=v2edcentcurrlm_codehigh68), alpha=0.5, color=NA) +
  scale_x_continuous(breaks=seq(1945,2015,10), limits=c(1945,2021), expand=c(0,0)) +
  scale_y_continuous(breaks=seq(-4,2,1), limits=c(-4.25,2.25)) +
  guides(color=guide_legend(override.aes=list(fill=NA)),
         fill=guide_legend(nrow=1)) +
  theme(legend.position="bottom", legend.title=element_blank(),
        axis.text.x=element_text(size=rel(1.5)),
        axis.text.y=element_text(size=rel(1.5)),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        legend.text=element_text(size=rel(1.5)),
        strip.background=element_rect(fill="white"),
        strip.text.x=element_text(size=rel(1.5), color="black"),
        plot.title=element_text(size=rel(1.5), hjust=0.5, color="black"))


ggsave(figf_9, file="Figures/Figure F-9.pdf", height=6, width=8)



###########################################################
## Figure F-10: Distribution of coder confidence ratings ##
###########################################################
figf_10 <- vindoc_coder %>%
  select(contains("conf")) %>%
  pivot_longer(cols=everything(),
               names_to="var",
               values_to="value") %>%
  ggplot(., aes(x=value)) +
  theme(plot.title=element_text(size=12, hjust=0.5),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.title=element_blank(),
        legend.position="none",
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour = "black", fill=NA, size=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey"),
        panel.grid.minor=element_blank(),
        plot.caption=element_text(hjust=0, size=12)) +
  geom_boxplot() +
  scale_x_continuous(breaks=seq(0,1,0.1), limits=c(0,1))

ggsave(figf_10, file="Figures/Figure F-10.pdf", height=2, width=5)



######################################################################
## Figure F-11: Mean coder-year confidence ratings across variables ##
######################################################################
coder_lgdppc <- coder_event_out("lgdppc", 0, NA, "Log GDPpc")

coder_reg_trans <- coder_event_out("reg_trans", 1, c("Autocratic\nTransition", "No\nTransition", "Democratic\nTransition"), "Regime Transition")
coder_v2svdomaut_ord <- coder_event_out("v2svdomaut_ord", 1, c("Not\nAutonomous", "Semi\nAutonomous", "Autonomous"), "Domestic Autonomy")
coder_v2elreggov <- coder_event_out("v2elreggov", 1, c("No", "Yes"), "Regional Gov't")

coder_is_war <- coder_event_out("is_war", 1, c("No War", "War"), "Inter-state War")
coder_cv_war <- coder_event_out("cv_war", 1, c("No War", "War"), "Intra-state War")
coder_coup <- coder_event_out("coup", 1, c("No", "Yes"), "Successful Coup")


figf_11 <- (coder_reg_trans | coder_v2svdomaut_ord | coder_v2elreggov) /
  (coder_is_war | coder_cv_war | coder_coup) &
  ylab(NULL) &
  theme(plot.margin = margin(5.5, 5.5, 5.5, 6.5))

figf_11 <- wrap_elements(figf_11) +
  labs(tag = "Mean Coder Confidence") +
  theme(plot.tag = element_text(size=12, angle = 90),
        plot.tag.position="left")

ggsave(figf_11, file="Figures/Figure F-11.pdf", height=8, width=10)



#####################################################################
## Figure F-12: Mean indicator-year confidence ratings across time ##
#####################################################################
coder_year <- vindoc_coder %>%
  select(year, contains("_conf")) %>%
  group_by(year) %>%
  summarize(across(where(is.numeric), mean, na.rm=TRUE)) %>%
  pivot_longer(cols=-year, 
               names_to="var", 
               values_to="value") %>%
  mutate(var=str_replace(var, "_conf", "")) %>%
  filter(var!="mean") %>%
  left_join(., read.csv("Data/vindoc_indicator_labels.csv")) %>%
  filter(!grepl("\\(", var_lab)) %>%
  group_by(var) %>%
  mutate(mean_conf=mean(value, na.rm=T)) %>%
  arrange(desc(mean_conf)) %>%
  mutate(var=factor(var, levels=unique(.[[2]]), labels=unique(.[[4]])))

figf_12 <- ggplot(coder_year, aes(x=year, y=var, fill=value)) +
  theme(plot.title=element_blank(),
        plot.margin=unit(c(0.75, 0.75, 0.75, 0.75), "cm"),
        axis.title.y=element_blank(),
        axis.text.y=element_text(size=12),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12),
        legend.title=element_blank(),
        legend.key=element_rect(colour=NA, fill=NA),
        legend.text=element_text(size=12),
        legend.position="bottom",
        legend.key.width=unit(1, "cm"),
        panel.border=element_rect(colour="black", fill=NA, linewidth=1),
        panel.background=element_rect(fill="white"),
        panel.grid=element_line(color="grey")) +
  geom_tile() +
  scale_fill_gradientn(colors=c("blue4","green1","green4"),
                       limits=c(0.75,1)) +
  scale_x_continuous(breaks=seq(1945, 2020, 10),
                     limits=c(1945, 2021),
                     expand=c(0,0)) +
  scale_y_discrete(limits=rev)

ggsave(figf_12, file="Figures/Figure F-12.pdf", height=8, width=10)



###############################################################
## Figure G-13: Indoctrination potential in the media (2021) ##
###############################################################
figg_13 <- world_map %>% 
  ggplot(aes(fill=v2xedvd_me_inco, map_id=region, alpha=is.na(v2xedvd_me_inco))) +
  geom_map(map=world_map, col=NA) +
  coord_fixed(1.3) +
  expand_limits(x=world_map$long, y=world_map$lat) +
  ggthemes::theme_map() +
  scale_fill_distiller(palette="Reds", guide="colorbar", direction=1,
                       na.value="gray60",
                       limits=c(0,1),
                       breaks=c(0,0.5,1)
  ) +
  scale_alpha_manual(values = c("TRUE" = .5), labels = c("N/A")) +
  scale_y_continuous(limits=c(-50,90)) +
  theme(legend.text=element_text(size=rel(2.5)),
        legend.title=element_blank(),
        legend.position=c(0.05,0.015),
        plot.title=element_text(size=rel(2), hjust=0.5),
        plot.background=element_rect(color="black", fill="white"),
        plot.caption=element_text(hjust=0, size=rel(1.45)))

ggsave(figg_13, file="Figures/Figure G-13.pdf", height=6, width=12)



#######################################################################
## Figure G-14: Indoctrination potential in the media across regimes ##
#######################################################################
figg_13 <- regime_plot(c("v2xedvd_me_inco"),
               c("Indoctrination potential"),
               "vdem",
               0,1,0,1,0.25,1)

ggsave("Figures/Figure G-14.pdf", height=4, width=4)



###########################################################
## Figure G-15: Indoctrination potential in media (2021) ##
###########################################################
figg_15 <- coef_all_func("v2xedvd_me_inco")
ggsave("Figures/Figure G-15.pdf", height=12, width=10)



###############################################################
## Figure H-16: Indoctrination potential in education (2021) ##
###############################################################
figh_16 <- coef_all_func("v2xed_ed_inpt")
ggsave("Figures/Figure H-16.pdf", height=12, width=10)



########################################################################
## Figure H-17: Democratic indoctrination content in education (2021) ##
########################################################################
figh_17 <- coef_all_func("v2xed_ed_dmcon")
ggsave("Figures/Figure H-17.pdf", height=12, width=10)



#######################################################################
## Figure H-18: Patriotic indoctrination content in education (2021) ##
#######################################################################
figh_18 <- coef_all_func("v2xed_ed_ptcon")
ggsave("Figures/Figure H-18.pdf", height=12, width=10)



###########################################################################
## Figure I-19: Levels of democracy and indoctrination potential in 2021 ##
###########################################################################
figi_19 <- dem_splot("v2xed_ed_inpt", "Indoctrination Potential Index")
ggsave("Figures/Figure I-19.pdf", height=9, width=9)



#####################################################################
## Figure I-20: Levels of democracy and democratic content in 2021 ##
#####################################################################
figi_20 <- dem_splot("v2xed_ed_dmcon", "Democratic Content Index")
ggsave("Figures/Figure I-20.pdf", height=9, width=9)



####################################################################
## Figure I-21: Levels of democracy and petriotic content in 2021 ##
####################################################################
figi_21 <- dem_splot("v2xed_ed_ptcon", "Patriotic Content Index")
ggsave("Figures/Figure I-21.pdf", height=9, width=9)



###################################################################
## Figure I-22: Democracy and the indoctrination indices in 2010 ##
###################################################################
figi_22 <- vindoc %>%
  mutate(e_p_polity=ifelse(-10>e_p_polity, 2, e_p_polity)) %>%
  select(country_id, year, v2x_liberal, v2x_polyarchy, e_p_polity, e_uds_median, v2xed_ed_inpt, v2xed_ed_dmcon, v2xed_ed_ptcon) %>%
  filter(year==2010) %>%
  select(-c(country_id, year)) %>%
  rename("Liberal Democracy"=1,
         "Electoral Democracy"=2,
         "Polity 5"=3,
         "UDS"=4,
         "Indoctrination Potential"=5,
         "Democratic Content"=6,
         "Patriotic Content"=7) %>%
  gather(key="dem_var", value="dem_value", 1:4) %>%
  gather(key="vindoc_var", value="vindoc_value", 1:3) %>%
  mutate(dem_var=factor(dem_var, levels=c("Liberal Democracy",
                                          "Electoral Democracy",
                                          "Polity 5",
                                          "UDS")),
         vindoc_var=factor(vindoc_var, levels=c("Indoctrination Potential",
                                                "Democratic Content",
                                                "Patriotic Content"))) %>%
  ggplot(., aes(x=dem_value, y=vindoc_value)) +
  theme(panel.background=element_blank(), 
        axis.line=element_line(color="black"),
        plot.title=element_text(size=12, hjust=0.5, color="black"),
        axis.text=element_text(size=12, color="black"),
        axis.text.y=element_text(size=12, color="black"),
        strip.text=element_text(size=12, color="black"),
        axis.title.y=element_blank(), 
        axis.title.x=element_blank(), 
        legend.title=element_blank(), 
        legend.key=element_blank(), 
        legend.position="none",
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.spacing=unit(0.75, "cm")) +
  geom_point() +
  geom_smooth() +
  stat_cor(aes(label=..r.label..)) +
  scale_y_continuous(limits=c(0,1.1), breaks=c(0,1)) +
  facet_grid(vindoc_var~dem_var, 
             scales="free")

ggsave(figi_22, file="Figures/Figure I-22.pdf", height=8, width=10)



###############################################################################################################
## Figure I-23: Indoctrination potential and content in education across regimes (BMR regime classification) ##  
###############################################################################################################
figi_23 <- regime_plot(c("v2xed_ed_inpt",
                   "v2xed_ed_dmcon", 
                   "v2xed_ed_ptcon"),
                 c("Indoctrination potential",
                   "Democratic indoctrination content",
                   "Patriotic indoctrination content"),
                 "bmr",
                 0,1,0,1,0.25,3)

ggsave("Figures/Figure I-23.pdf", height=4, width=12)



###################################################################
## Figure J-24: Components of the indoctrination potential index ##
###################################################################
figj_24 <- regime_plot(c("v2xed_ed_inco", 
                   "v2xed_ed_poed"),
                 c("Indoctrination coherence in education",
                   "Political values and ideology in education"),
                 "vdem",
                 0,1,0,1,0.25,2)

ggsave("Figures/Figure J-24.pdf", height=4, width=8)



############################################################################
## Figure J-25: Components of the democratic indoctrination content index ##
############################################################################
figj_25 <- regime_plot(c("v2edcritical", 
                   "v2edideolch_rec", 
                   "v2edplural", 
                   "v2edpoledrights"),
                 c("Critical engagement with education content",
                   "Ideology character in the curriculum",
                   "Pluralism in the curriculum",
                   "Political rights and duties in the curriculum"),
                 "vdem",
                 -2,2,-2.25,2.25,1,2)

ggsave("Figures/Figure J-25.pdf", height=8, width=8)



###########################################################################
## Figure J-26: Components of the patriotic indoctrination content index ##
###########################################################################
figj_26 <- regime_plot(c("v2edpatriot", 
                   "v2edscpatriot"),
                 c("Patriotic education in the curriculum",
                   "Presence of patriotic symbols in schools"),
                 "vdem",
                 -1.5,1.5,-1.5,1.6,1,2)

ggsave("Figures/Figure J-26.pdf", height=4, width=8)



###################################################################################
## Figure J-27: Components of the political education efforts in education index ##
###################################################################################
figj_27 <- regime_plot(c("v2edideol", 
                   "v2edpoledprim",
                   "v2edpoledsec"),
                 c("Ideology in the curriculum",
                   "Political education, primary school",
                   "Political education, secondary school"),
                 "vdem",
                 -1.5,1.5,-1.5,1.5,1,2)

ggsave("Figures/Figure J-27.pdf", height=8, width=8)



################################################################################
## Figure J-28: Components of the indoctrination coherence in education index ##
################################################################################
figj_28 <- regime_plot(c("v2xed_ed_cent", 
                   "v2xed_ed_ctag"),
                 c("Centralization of curriculum and textbooks",
                   "Control over educational agents"),
                 "vdem",
                 0,1,0,1,0.25,2)

ggsave("Figures/Figure J-28.pdf", height=4, width=8)



#####################################################################################
## Figure J-29: Components of the centralization of curriculum and textbooks index ##
#####################################################################################
figj_29 <- regime_plot(c("v2edcentcurrlm", 
                   "v2edcenttxbooks"),
                 c("Centralized curriculum",
                   "Centralized textbook approval"),
                 "vdem",
                 -2,2,-2.25,2.25,1,2)

ggsave("Figures/Figure J-29.pdf", height=4, width=8)



#####################################################################################
## Figure J-30: Components of the centralization of curriculum and textbooks index ##
#####################################################################################
figj_30 <- regime_plot(c("v2edteautonomy", 
                   "v2edteunionindp", 
                   "v2edtehire", 
                   "v2edtefire"),
                 c("Teacher autonomy in the classroom",
                   "Independent teacher unions",
                   "Teacher hiring",
                   "Teacher firing"),
                 "vdem",
                 -3,3,-3,3,1,2)

ggsave("Figures/Figure J-30.pdf", height=8, width=8)



#########################################################
## Figure K-31: Composite indoctrination index in 2021 ##
#########################################################
figk_31 <- world_map %>% 
  ggplot(aes(fill=v2xed_ed_indcomp, map_id=region, alpha=is.na(v2xed_ed_indcomp))) +
  geom_map(map=world_map, col=NA) +
  coord_fixed(1.3) +
  expand_limits(x=world_map$long, y=world_map$lat) +
  ggthemes::theme_map() +
  scale_fill_gradient2("Potential: ", 
                       low="darkred",
                       mid="yellow",
                       high="darkgreen",
                       guide="colorbar", 
                       na.value="gray60",
                       limits=c(-0.5,0.5),
                       breaks=c(-0.5,0,0.5)
  ) +
  scale_alpha_manual(values = c("TRUE" = .5), labels = c("N/A")) +
  scale_y_continuous(limits=c(-50,90)) +
  theme(legend.text=element_text(size=rel(2.5)),
        legend.title=element_blank(),
        legend.position=c(0.05,0.015),
        plot.title=element_text(size=rel(2), hjust=0.5),
        plot.background=element_rect(color="black", fill="white"),
        plot.caption=element_text(hjust=0, size=rel(1.45)))

ggsave(figk_31, file="Figures/Figure K-31.pdf", height=6, width=12)



#################################################################
## Figure K-32: Levels of democracy and indoctrination in 2021 ##
#################################################################
plot_dat <- vindoc %>%
  select(country_text_id, year, v2x_libdem, v2xed_ed_indcomp) %>%
  filter(year==2021) %>%
  select(-year) %>%
  rename("Liberal Democracy Index"=v2x_libdem,
         "Indoctrination Index"=v2xed_ed_indcomp)

figk_32 <- ggplot(plot_dat, aes(x=`Liberal Democracy Index`, y=`Indoctrination Index`, label=country_text_id)) +
  geom_hline(yintercept=0, linetype="dashed") +
  geom_point() +
  geom_text_repel(box.padding=0.2, point.padding=0.1, size=4, max.overlaps=30) +
  scale_y_continuous(limits=c(-0.5,0.5), breaks=seq(-0.5,0.5,0.5)) +
  theme_bw() +
  theme(axis.title.x=element_text(size=14),
        axis.title.y=element_text(size=14),
        axis.text.x=element_text(size=14, color="black"),
        axis.text.y=element_text(size=14, color="black"),
        plot.margin=unit(c(0.15,7,0.1,0.1), "cm"))  +
  annotate("text", x=0.045, y=0.5, label=paste0("Corr=", round(cor(plot_dat[,2], plot_dat[,3], use="pairwise.complete.obs"),2)), size=6) +
  coord_cartesian(xlim=c(0,1), clip="off") +
  annotation_custom(linesGrob(arrow=arrow(type="closed",
                                          ends="both",
                                          length=unit(3,"mm"))), xmin=1.1, xmax=1.1, ymin=-0.55, ymax=0.55) +
  annotate("text", x=1.25, y=0.475, label="Strong Democratic \n Indoctrination \n (High Democratic Content \n and Potential)", size=4.5) +
  annotate("text", x=1.25, y=0, label="Weak Indoctrination \n (Low Ideological Content \n and/or Potential)", size=4.5) +
  annotate("text", x=1.25, y=-0.475, label="Strong Autocratic \n Indoctrination \n (High Autocratic Content \n and Potential)", size=4.5)

ggsave(figk_32, file="Figures/Figure K-32.pdf", height=8, width=12)



##################################################
## Appendix L: Correlations and classifications ##
##################################################
## Load validation variables
validation_cy <- readRDS("Data/validation_cy.rds")
validation_ty <- readRDS("Data/validation_ty.rds")

## Get list of matches
validation <- read.csv("Data/validation_in.csv")

## Compute match statistics for each pair of matches
for (i in 1:nrow(validation)){
  if (validation$cy[i]==1){
    match_dat <- validation_cy
  } else{
    match_dat <- validation_ty
  }
  
  if (!is.na(validation$corr[i])){
    match_dat <- match_dat %>%
      select(country_id, year, all_of(validation$vindoc_var[i]), all_of(validation$match_var[i]), all_of(paste0(validation$vindoc_var[i],"_nr"))) %>%
      na.omit %>%
      rename(coder_ct=5)
    
    match_dat_3coder <- match_dat %>%
      filter(coder_ct>=3)
    
    validation$corr[i] <- round(cor(match_dat[,3], match_dat[,4]),2)
    validation$corr_3coder[i] <- round(cor(match_dat_3coder[,3], match_dat_3coder[,4]),2)
    validation$country[i] <- length(unique(match_dat$country_id))
    validation$years[i] <- paste0(min(match_dat$year),"-",max(match_dat$year))
    validation$n[i] <- nrow(match_dat)
  } else{
    match_dat <- match_dat %>%
      select(country_id, year, all_of(paste0(validation$vindoc_var[i],"_ord")), all_of(validation$match_var[i]), all_of(paste0(validation$vindoc_var[i],"_nr"))) %>%
      na.omit %>%
      rename(coder_ct=5)
    
    match_dat_3coder <- match_dat %>%
      filter(coder_ct>=3)
    
    tab_dat <- table(match_dat[,3], match_dat[,4])
    tab_dat <- table(match_dat_3coder[,3], match_dat_3coder[,4])
    
    validation$class[i] <-  round(sum(diag(tab_dat))/sum(tab_dat),2)
    validation$class_3coder[i] <-  round(sum(diag(tab_dat))/sum(tab_dat),2)
    validation$country[i] <- length(unique(match_dat$country_id))
    validation$years[i] <- paste0(min(match_dat$year),"-",max(match_dat$year))
    validation$n[i] <- nrow(match_dat)
  }
}

summary(c(abs(validation$corr), abs(validation$class)))
summary(c(abs(validation$corr_3coder), abs(validation$class_3coder)))

validation_rank <- validation %>%
  select(vindoc_var, match_var, corr, class) %>%
  mutate(match=ifelse(!is.na(corr), abs(corr), abs(class))) %>%
  arrange(desc(match))

## Table 1
print(validation_rank)



###############################
## Appendix M: Linz analysis ##
###############################
## Table M-3
linz_summary_all <- linz_dat %>%
  select(v2xed_ed_inpt, v2xed_ed_inpt_mr,
         gwf_party, gwf_personal, gwf_military, gwf_monarch, lgdppc) %>%
  filter(!is.na(v2xed_ed_inpt)) %>%
  describe() %>%
  select(n, mean, sd, min, max) %>%
  as.data.frame
round(linz_summary_all, 4)

## Table M-4
linz_summary_3coder <- linz_dat %>%
  select(v2xed_ed_inpt, v2xed_ed_inpt_mr,
         gwf_party, gwf_personal, gwf_military, gwf_monarch, lgdppc) %>%
  filter(!is.na(v2xed_ed_inpt)) %>%
  filter(v2xed_ed_inpt_mr>=3) %>%
  describe() %>%
  select(n, mean, sd, min, max) %>%
  as.data.frame
round(linz_summary_3coder, 4)

## Table M-5
linz_ed_inpt_twoway <- linz_out(
  "v2xed_ed_inpt",
  "v2xed_ed_inpt_mr",
  "twoways")

stargazer(fig13$fe_reg, 
          fig13$fe_reg_3coder,
          linz_ed_inpt_twoway$fe_reg, 
          linz_ed_inpt_twoway$fe_reg_3coder, 
          type="text",
          dep.var.caption="DV: Indoctrination Potential in Education",
          dep.var.labels.include=F,
          column.labels=rep(c("All observations","At least 3 coders"), 2),
          covariate.labels=c("Party", "Personalist", "Monarchy", "GDPpc (log)"),
          keep.stat=c("n"),
          add.lines=list(c("Countries", 
                           rep(c(
                             length(unique(index(fig13$fe_reg, "country_id"))),
                             length(unique(index(fig13$fe_reg_3coder, "country_id")))), 
                             2)),
                         c("Country FE", 
                           rep("Yes", 4)),
                         c("Year FE",
                           rep(c("No", "Yes"), each=2))),
          no.space=T,
          notes.align="l",
          table.layout ="-lc#-t-sa-n",
          digits=4)

## Figure M-33 and M-34
linz_ed_inpt_codelow68 <- linz_out(
  "v2xed_ed_inpt_codelow68",
  "v2xed_ed_inpt_mr",
  "individual")

linz_ed_inpt_codehigh68 <- linz_out(
  "v2xed_ed_inpt_codehigh68",
  "v2xed_ed_inpt_mr",
  "individual")

linz_ed_poed <- linz_out(
  "v2xed_ed_poed",
  "v2xed_ed_poed_mr",
  "individual")

linz_ed_inco <- linz_out(
  "v2xed_ed_inco",
  "v2xed_ed_inco_mr",
  "individual")

linz_me_inco <- linz_out(
  "v2xedvd_me_inco",
  "v2xedvd_me_inco_mr",
  "individual")

linz_supp_all <- rbind(linz_ed_inpt_codelow68$plot_dat,
                       linz_ed_inpt_codehigh68$plot_dat,
                       linz_ed_poed$plot_dat,
                       linz_ed_inco$plot_dat,
                       linz_me_inco$plot_dat) %>%
  filter(subset=="All observations") %>%
  mutate(dv=rep(c("linz_ed_inpt_codelow68",
                  "linz_ed_inpt_codehigh68",
                  "linz_ed_poed",
                  "linz_ed_inco",
                  "linz_me_inco"), each=4)) %>%
  mutate(dv=factor(dv, 
                   levels=c("linz_ed_inpt_codelow68",
                            "linz_ed_inpt_codehigh68",
                            "linz_ed_poed",
                            "linz_ed_inco",
                            "linz_me_inco"),
                   labels=c("Indoctrination Potential in Education 68% CI Lower Bound",
                            "Indoctrination Potential in Education 68% CI Upper Bound",
                            "Political Education Effort in Education",
                            "Indoctrination Coherence in Education",
                            "Indoctrination Potential in the Media")))

linz_supp_3coder <- rbind(linz_ed_inpt_codelow68$plot_dat,
                          linz_ed_inpt_codehigh68$plot_dat,
                          linz_ed_poed$plot_dat,
                          linz_ed_inco$plot_dat,
                          linz_me_inco$plot_dat) %>%
  filter(subset=="Observations with at least 3 coders") %>%
  mutate(dv=rep(c("linz_ed_inpt_codelow68",
                  "linz_ed_inpt_codehigh68",
                  "linz_ed_poed",
                  "linz_ed_inco",
                  "linz_me_inco"), each=4)) %>%
  mutate(dv=factor(dv, 
                   levels=c("linz_ed_inpt_codelow68",
                            "linz_ed_inpt_codehigh68",
                            "linz_ed_poed",
                            "linz_ed_inco",
                            "linz_me_inco"),
                   labels=c("Indoctrination Potential in Education 68% CI Lower Bound",
                            "Indoctrination Potential in Education 68% CI Upper Bound",
                            "Political Education Effort in Education",
                            "Indoctrination Coherence in Education",
                            "Indoctrination Potential in the Media")))

figm_33 <- ggplot(linz_supp_all, aes(x=vars, y=coef, group=fct_rev(dv), color=dv)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey")) +
  geom_errorbar(aes(ymin=lb, ymax=ub, width=0), size=1, position=position_dodge(width=0.4)) +
  geom_point(position=position_dodge(width=0.4), size=2) +
  geom_hline(yintercept=0, linetype="dashed") +
  guides(color=guide_legend(direction="vertical")) +
  scale_y_continuous(limits=c(-0.2,0.14)) +
  coord_flip()

ggsave(figm_33, file="Figures/Figure M-33.pdf", height=6, width=6)


figm_34 <- ggplot(linz_supp_3coder, aes(x=vars, y=coef, group=fct_rev(dv), color=dv)) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_text(size=12, color="black"),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_text(size=12, color="black"),
        legend.position="bottom",
        legend.title=element_blank(),
        legend.key=element_rect(fill=NA),
        legend.text=element_text(size=12),
        panel.border=element_rect(colour="black", fill=NA, size=1),
        panel.background=element_rect(fill="transparent"),
        panel.grid.major=element_line(color="lightgrey"),
        panel.grid.minor=element_line(color="lightgrey")) +
  geom_errorbar(aes(ymin=lb, ymax=ub, width=0), size=1, position=position_dodge(width=0.4)) +
  geom_point(position=position_dodge(width=0.4), size=2) +
  geom_hline(yintercept=0, linetype="dashed") +
  guides(color=guide_legend(direction="vertical")) +
  scale_y_continuous(limits=c(-0.2,0.14)) +
  coord_flip()

ggsave(figm_34, file="Figures/Figure M-34.pdf", height=6, width=6)



###############################################################
## Figure N-35: Patriotism and nationalism in the curriculum ##
###############################################################
fign_35 <- vindoc %>% 
  ggplot(aes(x=v2xed_ed_ptcon, y=v2edideolch_1)) +
  geom_point(size=1) +
  labs(x = "Patriotism in language curriculum", 
       y = "Nationalist ideology in history curriculum",
       title = "Comparison: all years") +
  geom_smooth() +
  stat_cor(method = "pearson", label.y = .95, label.x = .01) +
  theme(panel.background=element_blank(),
        axis.line=element_line(color="black"),
        plot.title=element_text(size=rel(1), hjust=0.5, color="black"),
        axis.text=element_text(size=rel(1), color="black"),
  ) +
  
  vindoc %>% filter(year==2021) %>%
  ggplot(aes(x=v2xed_ed_ptcon, y=v2edideolch_1, label=country_text_id)) +
  geom_point(size=1) +
  labs(x = "Patriotism, language curriculum", 
       y = "Nationalist ideology, history curriculum",
       title = "Comparison: 2021") +
  stat_cor(method = "pearson", label.y = 0.9, label.x = .01) +
  geom_text_repel(size=2.5, max.overlaps=30, color="black") +
  theme(panel.background=element_blank(),
        axis.line=element_line(color="black"),
        plot.title=element_text(size=rel(1), hjust=0.5, color="black"),
        axis.text=element_text(size=rel(1), color="black"),
        axis.text.y=element_text(vjust=0.35),
  ) +
  scale_y_continuous(limits=c(-0.1, 1), breaks=seq(0, 1, 0.25)) +
  
  plot_layout(ncol=1)

ggsave(fign_35, file="Figures/Figure N-35.pdf", height=10, width=6)


## Close log file
log_close()
